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We demonstrate control over heat flow in an A''-terminal molecular junction. Using simple model 
Hamiltonians we show that the heat current through two terminals can be tuned, switched, and 
amplifled, by the temperature and coupling parameters of external gating reservoirs. We discuss 
two models: A fully harmonic system, and a model incorporating anharmonic interactions. For 
both models the control reservoirs induce thermal fluctuations of the transition elements between 
molecular vibrational states. We find that a fully harmonic model does not show any controllability, 
while for an anharmonic system the conduction properties of the junction strongly depend on the 
parameters of the gates. Realizations of the model system within nanodevices and macromolecules 
are discussed. 
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Control over vibrational energy flow in nanoscale struc- 
tures and single molecules is a long standing goal in 
many parts of physical science and nanotechnology. His- 
torically, intramolecular vibrational redistribution (IVR) 
vi^as a topic of great interest in chemistry and physics. 
IVR processes must be reckoned for understanding, and 
ultimately controlling, molecular dynamics and chemi- 
cal kinetics The efficiency of these processes is the 
basic assumption behind the well validated RRKM re- 
action rate theory 0, From a different perspective, 
the unexpected results of the computer experiment of 
Fermi-Pasta- Ulam [3], showing no equipartition of en- 
ergy among normal modes in harmonic chains including 
small nonlinear terms, lead to extensive research of IVR 
in nonlinear systems Q- 

Recurrent theoretical interest in this field is due to 
the impressive progress in probing thermal properties of 
nanoscale systems such as nanotubes 6, 7, Sj^self assem- 
bled monolayers 0, EH , and thin films [T^, and due 
to the development of more tunable systems [iJl • Recent 
progress in molecular electronics and nanomechanics has 
raised further interest in exploring mechanisms of energy 
flow in nano-level systems. In molecular electronics, local 
heating of nanoscale devices might cause structural in- 
stabilities undermining the junction integrity p^[l5l.[l^. 
Engineering good thermal contacts and cooling of the the 
junction are necessary for a stable operation mode. Min- 
imization of mechanical devices, e.g. refrigerators [TtI 
and pumps [31 , to the molecular scale is a topic of great 
interest for technologies such as chemical sensing, power 
generation and energy conversion [H, [H, [13, [IH . In this 
context, it is crucial to understand, and ultimately con- 
trol, the dynamics of phonons in nanoscale structures, or 
analogously, vibrational modes in molecular systems. 

The heat conduction properties of molecular chains can 
be tuned cither by manipulating the internal molecular 
parameters, or by externally gating the system. We re- 
fer to the flrst mode of control as " static" , or " internal" , 
while, as we explain below, we consider gating as a "dy- 
namical", or an "external" control scheme. 




FIG. 1: A schematic representation of the model studied in 
this work, exemplified on a three-bath system. The molecular 
units (central part) interacts with multiple thermal reservoirs 
maintained at constant temperatures. Heat current is defined 
as positive when flowing into the central molecular segment. 



Static control is realized by adjusting internal system 
parameters, e.g. atomic masses and interatomic potential 
energy, with the goal of increase/decrease of the system 



thermal conductivity [24I . This type of control problem 



has been extensively discussed in the context of Fourier 
transport. Here the main challenge is to identify the 
necessary and sufficient conditions for the validity of the 
Fourier law of heat conductivity, J cx —"VT, in low di- 
mensions and for quantum systems [H, [l^, [2^, [2^ . By 
Engineering the molecular system, one can also build 
functional devices, for example, a thermal rectifier, where 
the nanojunction allows heat flux in one direction, while 
it acts like an insulator when the temperature gradient 
is reversed [13, 113, [H, [H ■ 

In this context we would like to emphasize that tuning 
the thermal conductance of harmonic chains is also fea- 
sible, though these systems demonstrate abnormal (non- 
Fourier) thermal current. In the language of the thermal 
Landauer formula (30l |. this can be accomplished by ma- 
nipulating the transmission coefficient for phononic heat 
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flow through the device, e.g. by introducing impurities 



into the structure [31 1 



In this paper we present a simple model that illustrates 
an external control over thermal current in molecular sys- 
tems. The generic setup includes a molecule (subsystem) 
coupled to N thermal baths of difi^erent temperatures. 
Heat current flowing through the system may be modi- 
fied by a control reservoir, where in general, more than 
one gate may couple to the subsystem. Typical control 
parameters are the gate temperature and its coupling 
to the molecular unit. A schematic representation of a 
three-bath scenario is shown in Fig. [T] 

We consider two realizations of this standard setup: (i) 
Artificial nanodevices operating as thermal amplifiers or 
transistors [3^ . Here the generic system includes three 
segments, source (L) drain (i?) and gate (G), following 
the notation used in semiconductor transistors, where 
heat flow between the source and drain is controlled by 
the gate terminal. This device may be realized by fabri- 
cating branched nanotubes with T, Y, and X shapes [ssj . 
(ii) Macromolecules, e.g. proteins and dendriemds with 
spatially separated sidegroups. Energy may be funneled 
between the molecular groups by a control unit (G), for 
example a protein, that can temporarily bind to the sys- 
tem. We refer to the energetically hot group as L, while 
the R group is the energy accepting sidegroup. 

The role of the gate (G) may be modeled as inducing 
thermal fluctuations of the L - R transition element [s^ , 
thus we refer to this mode of control as dynamical. The 
control element can also be identified as the solvent it- 
self, modulating system parameters. Inserting a molecule 
in different mediums may therefore modify its heat con- 
duction properties. In this context this work provides a 
simple framework for investigating IVR in solutions [35| . 

The main question to be addressed in this paper is 
what is the role of anharmonic interactions, specifically, 
nonbilinear molecule-surface interactions, in controlling 
the thermal properties of a gated system. It is widely 
accepted that nonlinear interactions are essential for 
sho wing normal (Fourier) transport in molecular chains 
[2^ l36l |. Anharmonic interactions are also necessary 
for bringing in rectifying behavior [13, [H, [2^, and for 
manifesting nonlinear thermal conductance characteris- 
tics [32,, ,37,1 . We show next that anharmonic interactions 
are also a crucial element for realizing dynamical control 
of heat current. 

We consider two models. The first system is a proto- 
type for transport in harmonic chains. The second model 
incorporates anharmonic interactions in the molecule, 
and also assumes nonbilinear system-bath couplings. For 
both models we focus on two quantities: (i) We calculate 
the heat current at the terminal u, J^, and (ii) we inves- 
tigate the net heat current flowing between two surfaces, 
Ju^^i- The objective of our calculation is to demonstrate 
that these quantities strongly depend on the parameters 
of the gate reservoirs (temperature and energetics) for 
the anharmonic model only. 

We begin with the harmonic model. In this case 



both the reservoirs (inverse temperatures Pi, = ^ 
V ~ 1...7V) and the molecular unit are modeled by a set of 
non interacting bosons coupled via a bilinear term. For 
simplicity, we assume that heat transfer in the molecular 
unit is dominated by a speciflc single mode. 
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+ (6t + 6)^A,-.(a],, + a,-,). (1) 



Here 6^ (6) are creation (annihilation) operators for the 
molecular mode of frequency wq. Similarly, a j ^ (oj,!/) 
are the operators for the mode j of the v reservoir. Since 
system-bath interaction is bilinear, the model Hamilto- 
nian can be exactly diagonalized, to be represented in 
terms of a set of noninteracting phonons. We refer to 
model HI) as the boson-boson (b-b) Hamiltonian. 

The dynamics of the subsystem can be exactly solved 
using various techniques, e.g. the generalized Langevin 
equation (3ll . [sH] and Master equation formalism [H, 113, 
l38j . The result of these calculations is the " thermal Lan- 
dauer formula" [l^, [3l| , where in the classical limit the 
heat current (for N = 2) linearly depends on the tem- 
perature difference between the two thermal baths. We 
briefly follow here the derivation within the Master equa- 
tion formalism, generalizing the results of Refs. 29, 33] 
for an N terminal system. 

Under the assumption of weak system-bath interac- 
tions, going into the Markovian limit, the probabilities 
Pn to occupy the n state of the molecular oscillator sat- 
isfy the Master equation [53, [11] , 

Pn = nkuPn-i{t) 

+ {n + l)kdPn+i{t) - [nkd + {n + l)ku]Pnit).{2) 

The nonadiabatic relaxation and excitation rates, kd and 
ku respectively, are given by summing up contributions 
from each reservoir, as no correlations exist between the 
different baths {ly — 1..N), 



ku — / ^ ; kd — > k^ 



(3) 



It can be shown that 



kd ^T^{LUo){l + n^{ujo)); fcj^ = r^(LJo)"i/(t^o), (4) 
where 

T^{u!o) = 2Tr'^X'^jS{ujj - ujq). (5) 

Here n^{uj) = [e^/'^" — 1]"^ is the Bose- Einstein distri- 
bution function for the v reservoir. The heat current 
properties of this model are obtained from the steady 
state solution of Eq. ^ with the nonadiabatic rates 
. The steady state heat flux at the terminal is given 
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by calculating the difference between heat flow from the 
V bath into the the molecular mode, leading to vibra- 
tional excitations within the molecule, and the outgoing 
molecule- reservoir energy current, resulting from relax- 
ation processes inside the molecule. 



(6) 



The current is defined positive when flowing from the 
contact into the molecule. In the classical limit {Ty > loq) 
we get [13 



(7) 



Considering this expression, we can identify the directed 
current v ^ ii as 



j{bb) _ 



(8) 



We refer next to two specific terminals as source {L) and 
drain (i?), while all other N — 2 baths are referred to as 
gates (G). When currents in the G terminals are zero, 
Jv^L,B. ~ 0, the gates acquire the same temperature in 
steady state, Tq = [TlTl + TrTr)I{Tl + Tr). The 
current at the L/R contact is then given by 



J 



ibb) 



Li R 



r, 



(9) 



which is the same result as obtained when Tg = 0, see 
Eq. (O. We also find that the current Jl^r decays with 
the number of thermal reservoirs as 



J 



(6b) 
L^R 



TlTr 



{Tl^Tr); (Jg=0), 



(10) 



due to additional decay channels. To summarize, we find 
that in the harmonic limit the effect of gate terminals 
is simply to effectively increase the broadening F, while 
the gates' temperatures can not modify the current in the 
system. Thus, there is no control over the heat dissipated 
(or absorbed) from the contacts in a purely harmonic 
system. 

Next we show that in a model consisting nonlinear in- 
teractions heat current can be strongly controlled by the 
temperature of a gate terminal. As a case study we con- 
sider the spin-boson (s-b) model, generalized to include 
N bosonic reservoirs (creation operators Oj y, ly — 1...N) 
linearly coupled to a spin (two-level) system. 



(11) 
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Here Wq is the energy difference between the spin lev- 
els with tunneling splitting A. In this model internal 



molecular anharmonicity is introduced by truncating the 
spectrum of the single molecular mode to include only 
the lowest two energy states. We do not allow for other 
phonon-phonon scattering processes, e.g. umklapp pro- 
cesses, that can lead to normal conductivity as in the 
Peierls model [1^. Using the small polaron transforma- 
tion it can be shown that this model represents a 
molecular mode coupled nonbilinearly to the harmonic 
reservoirs (291, 
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Je'''a++e-'''a-)+HsMff (12) 

Here ^ J^^ ^j^; ^-^ = «Ej l^(al,i/ - o-3,i')^ Hsh^ft = 
ly 4^'." is an energy shift henceforth incorporated into 
the zero order energies. Eq. (fT2|) shows that the role of 
the thermal reservoirs is to modulate the transition ele- 
ments between molecular vibrational levels. The impor- 
tant feature of this model is that system-bath couplings 
(Eq. are multiplicative in the bath degrees of free- 

dom, rather than additive as is the linear harmonic model 
([3]) . We do not distinguish in this model between the role 
of the different reservoirs (source, drain and gates). One 
could construct variants of this model, where the gates 
interact in a distinct functional form. For small A the 
Hamiltonian leads again to nonadiabatic dynamics, Eq. 
0, with 71 = 0, 1. The and rate constants are given by 



A2 A2 

ku = —C{-ujo); kd = —C(ujq) 



Here 



C{t) 



(13) 



dte*"«*C(t), 
n,a(t); C,{t) = (e^f^.We-if^^W),. (M) 



The trace is performed over the v reservoir degrees of 
freedom. For convenience, in what follows we disregard 
the prefactor (A/2)^. Using the convolution theorem, it 
can be shown that the function C(u;) can be decomposed 
in terms of the N reservoirs correlation functions. 



C(^o) 



duj2 



du 



N-l 



xCl{LJi)C2{iO2)---CN{uJ0 — UJl — UJ2 



LUN-l) 

(15) 



where Cy{uj) — 



*Ci,{t)dt is identified as the rate 



constant affected from the v thermal bath. The heat 
flux at the v terminal can be formally written for an 
iV terminal system by considering combinations of all 
processes in which the reservoirs exchange energy with 
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the subsystem [29|, 

/oc 
uj^dujiduj2---duJi,...dujN-i^k^i^,NCkiuJk) x 
-oo 

Cu(-UJu)Cn{-(^1 - t^2--- + i^y- - i^N-l - Wo)^o] ■ 

(16) 

For clarity, we include the explicit expression for the heat 
flux in a three-terminal junctions, measured at terminal 
'1', 

j['''\N = -i) = 

/OO 
dL02[Ci{(jJi)C2{i02)Cz{(J^i^ — iOi— L02)Pl 
-oo 

-Ci{-uJi)C2{oJ2)Cz{-uJa +0Ji- oJ2)Pa\- (17) 

The population of the spin levels is given by Pq = 1— Pi — 
C(a;o)/[C(tJo) + C'(— Wo)]- Assuming strong couphng, go- 
ing into the high temperature classical limit, (T^ > ujq), 
the kernel C{t) can be calculated in the short time limit 

a(i) = e-*-W; <P,{t)^iEl,t + T,Elit\ (18) 

where E'^j — is the reorganization energy of the 

u reservoir. In frequency domain we find that 

{uj^-Em?^ 



expression, leading to strong controllability, as opposed 
to the harmonic results, Eqs. ([7]) and ([5]). 

We exemplify control over the heat current in the sys- 
tem by studying a source-drain-gate situation, where two 
reservoirs are considered as source {L) and drain (i?), 
while N — 2 baths are identified as gates (G) . Under the 
condition of zero current in the gating terminals ( Jc = 
0) , their temperatures can be determined self consistently 
to yield To - Tm = (^^f/Ti + E'^jTr) / {E'^, + Eg). The 
current at the source/drain contact then becomes 



Ciioo) - 



TT 

TmEm 

1 



exp 



^TmE 



: exp 



M 



^T.El, 



(19) 



where Em is the total (TV baths) reorganization energy 
and Tm is an effective temperature for the subsystem, 



T, 



Em - ^Em; 

J2u EmT„ 



M 



Em 



(20) 



Integrating Eq. p6)) utilizing (|T9|) yields the heat current 
at the v contact (ignoring a multiplicative numeric factor 
of %/47r). 



M 



(i?MTM)3/2(l + e-o/TM)' 



(21) 



The current flowing through the system between the ter- 
minals v and /i is given by 



'^'^ ^' e ^i^mTm . (22) 



'^^^'^^ (i?MTM)3/2(l + e-o/TM) 



The temperatures of the gating terminals and their cou- 
plings to the molecule appear in a nontrivial way in this 



j: 



(sb) 



{sb) 



-J 



in - tr 



Em Em 



E 



M 



+ Eg [EMTcf^ (1 + e-o/^o) ■ 



(23) 



When all gates evenly couple to the subsystem, E'^j = 
E^j, and for Em = NE^j > wq, we find that the current 
exponentially decays with TV (corrected by a power law) , 
jf^ oc AriV-i/2e-JVS°,/4TG^ AT = Tl- Tr. We can 
also calculate the L ^ R current, again taking = 0, 



J^l1r = {Tl-Tr) 



^ T\T ^ 



{EmTgY'^ (1 + e-o/To)- 



(24) 



In the limit of strong coupling, NE^j > ujq, we find 
that cx ATiV-3/2e-^^M/4TG, Therefore, the 

temperature Tq serves as an effective activation tem- 
perature, exponentially enhancing the directed current, 
while Em is the energy gap for transport. Note that 
Jl/Jl^r = Em /{Eg + Eg), in analogy with the be- 
havior of the fully harmonic model, Eq. ((9|)- (fT0|) . 

We present next numerical results calculated by ap- 
plying Eqs. (PT1) - (P^ on a three-terminal (L, R and G) 
configuration. The control parameters are the coupling 
strength of the gate to the subsystem, given in terms of 
reorganization energy Eg and the temperature of the G 
reservoir. For simplicity we take A/2 = 1 meV. 

Figure [2] displays results in the absence of (net) energy 
flow between the gate and the subsystem. Taking all 
coupling strengths to be equal sets the gate temperature 
to Tg = (Tl + Tr)/2. Simulating Eq. ^ we find 
that the current at the L/R interfaces strongly decays 
with Eg. We also show strong control over the directed 
L ^ R current using Eq. ([M]) . 

Next we modulate the temperature of the gate reser- 
voir, and manifest that the system can act like a good 
thermal conductor, or an insulator, depending on the 
gate temperature. Figure [3] shows that the system be- 
haves as a bad thermal conductor at low gating temper- 
ature Tg ~ lOOK, while at Tg = 250^ (where Jg = 0) 
the current is amplified by a factor of ~10. 

Finally, in Figure [J] we analyze the directed heat flux 
Jl^r as a function of the gate temperature. We show 
strong amplification of current, of three orders of mag- 
nitudes, when the gate terminal is strongly coupled to 
the molecular mode (dashed-dotted line). In contrast, at 
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FIG. 2: Control over heat current in a three-terminal config- 
uration under the condition of zero current at the gate Jg=0. 
The heat current at the L contact (full) strongly decays with 
E'^j. Also shown are the directed currents Jl^r (dashed), 
Jg^r (dashed-dotted), and Jg-^l (dotted). Other param- 
eters are Tl=300K, rR=200K, Tg=250K, EIi = E'^i = 50 
meV. 
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FIG. 3: Switching the heat current in a three-terminal setup: 
The system behaves as an insulator for low control temper- 
atures, while for high Tg the system becomes a good ther- 
mal conductor. Jl (dashed) and —Jr (full) are equal at 
Tg = 250if, where Jg (dotted) diminishes. Other parameters 
are Tl=300K, rR=200K, E^i = Sf^ = 50 meV, E% = 300 
meV. 

weak coupling (full line), Jl^r is insensitive to the pres- 
ence of the G channel. As discussed above, the model 
Hamiltonian pT|) could be also introduced for describing 
transport between two sidegroups in a macromolecule im- 
mersed in a solvent. In this context, Fig. 2] reveals that 
the solvent may enhance the L to R current, in analogy 
with the effect of solvent assisted IVR [s^ . 

In conclusion, we have presented here two models for 
the study of thermal transport in gated systems: a purely 



FIG. 4: Amplification of the heat current in a three-terminal 
configuration. The current Jl^r can be strongly manip- 
ulated by the control bath at strong coupling. Tl =300K, 
rH=200K, Eii = E^i = 50 meV, £^=10 meV (full), 100 
(dashed), 200 (dotted) and 300 meV (dashed-dotted). 



harmonic system, and an anharmonic model. In both 
cases the gate terminals induce thermal fluctuations of 
the transition elements between the molecular vibra- 
tional states, leading to dynamical control of heat cur- 
rent. We (trivially) found that the purely harmonic setup 
cannot bring in a gated behavior. In contrast, in the spin- 
boson model, incorporating anharmonic interactions, the 
system can behave either as an insulator, or as a good 
thermal conductor, depending on the gates parameters: 
temperatures and molecule-bath coupling strength (bind- 
ing energy using proteins terminology). 

We may also explore other variants of the anharmonic 
model. For example, the central molecular unit can 
bilinear ly- weakly couple to the L and R surfaces, while 
the gates may couple nonbilinearly-strongly to the sub- 
system. Such models should basically display the same 
characteristics as discussed above. 

The effects described in this paper may be also stud- 
ied u sing classical molecular dynamics (MD) simula- 
tions |32l |. The advantage of our formalism over such 
Langevin equation treatment is twofold: (i) Nonbilinear 
system-surface couplings are difficult to simulate within 
Langevin dynamics, leading to a coordinate dependent 
friction coefficient [4l|. (ii) The net heat current be- 
tween two reservoirs, J^^fi, cannot be resolved within 
MD simulations, only the current at each terminal can 
be directly obtained. In contrast, one can clearly identify 
these currents in the analytical expression (|22|) . 

The models presented in this paper can be realized by 
fabricating an N terminal nanodevice. The heat con- 
ductance of suspended two-terminal nanotubes was mea- 
sured by detecting changes in the electric resistance of 
the attached heater and detector pads [l^. We suggest 
employing this method in a three-terminal setup, e.g. 
by connecting a T-shaped nanotube to three conducting 
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surfaces. Macromolecules with spatially localized side- 
gro up^ also offer a playground for realizing these models. 

Nonlinear electrical devices (rectifiers, switches and 
transistors) have shaped technology in the last 60 years. 
Nonlinear nanomechanical devices promise to revolution- 
ize the technology of the future as well, whereas phonons. 



instead of electrons become the carriers of information 
and the computation element. This paper, presenting a 
simple study for the control of heat flow in nanosystems, 
is a first step in addressing this challenge. Introducing 
quantum effects will further offer new pathways for en- 
ergy control at the nanoscale [i^ . 
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